#include <stdio.h>
#include <malloc.h>
#include <float.h>
#include <math.h>
#include "clam.h"

#define MAX_PRECOMP 250  /*WMN 6/05 because it gets too slow to make the array */

extern void choose_sulston(struct tmpCz* clone1, struct tmpCz* clone2, struct tmpSz* result);

static double*** sulston_precomp = 0;
static int precomp_max_bands;
static int old_tolerance;

void precompute_sulston();
void destroy_precomp_array();
void evaluate_sulston(int b1, int b2);

/* Get the precomputed value, IF there is one. There can fail to be one if
precomputation was not turned on (in which case precomp_max_bands = 0),
or if the clone has too many bands (the only way I can think for this 
to happen is if the array is precomputed, and then new clones are added
to the project, which have higher band counts). */

int get_precomp_value(int b1, int b2, int m, double* result)
{
    if (precomp_max_bands == 0 || b1 > precomp_max_bands ||  b2 > precomp_max_bands) 
    {
        return 0;
    }  
    if (b1 >= b2)
    {
        *result = sulston_precomp[b1][b2][m];
    }
    else
    {
        *result = sulston_precomp[b2][b1][m];
    }
    return 1;
}
int get_max_band_count()
{
   int i, max = 0;
   struct tmpCz Cx;;
   for (i = 0; i < arrayMax(acedata); i++)
   {
      loadCzfast(&Cx,i);
      if (Cx.nbands > max)
      {
         max = Cx.nbands;
      }
   }
   return max;
}

void destroy_precomp_array()
{
    int b1, b2;
    if (sulston_precomp == 0)
    {
        return;
    }
    for (b1 = 1; b1 <= precomp_max_bands; b1++)
    {
        if (sulston_precomp[b1]== 0)
        {
            break;
        }
        for (b2 = 1; b2 <= b1; b2++)
        {
            if (sulston_precomp[b1][b2] == 0)
            {
                break;
            }
            free(sulston_precomp[b1][b2]);
            sulston_precomp[b1][b2] = 0;
        }
        free(sulston_precomp[b1]);
        sulston_precomp[b1] = 0;
    }
    free(sulston_precomp);
    sulston_precomp = 0;
    precomp_max_bands = 0;
}
void create_precomp_array(int max_bands)
{
    int b1, b2;
    int allocated = 0;

    precomp_max_bands = 0; /* so if malloc fails, precomp will not be used */
    sulston_precomp = (double***)malloc( (max_bands + 1)*sizeof(double**));
    if (sulston_precomp == 0)
    {
        printf("Unable to allocate memory for precomputation array\n");
        return;
    }
    allocated += (max_bands + 1)*sizeof(double**);
    for (b1 = 1; b1 <= max_bands; b1++)
    {
        sulston_precomp[b1] = (double**)malloc( (b1 + 1)*sizeof(double*));
        if (sulston_precomp[b1] == 0)
        {
            printf("Unable to allocate memory for precomputation array entry %d\n", b1);
            return;
        }
        allocated += (b1 + 1)*sizeof(double*);
        for (b2 = 1; b2 <= b1; b2++)
        {
            sulston_precomp[b1][b2] = (double*)malloc( (b2 + 1)*sizeof(double));
            if (sulston_precomp[b1][b2] == 0)
            {
                printf("Unable to allocate memory for precomputation array entry %d,%d\n", b1,b2);
                return;
            }
            allocated += (b2 + 1)*sizeof(double);
        }
    }
    precomp_max_bands = max_bands;
    allocated /= 1024;
    printf("sulston precomputation: allocated %d kb\n",allocated);
}
void precompute_sulston()
{
    int b1, b2, m;
    int b1_save, b2_save, m_save, max_bands;

    if (0 == Pz.precompute)
    {
        return;
    }

    max_bands = get_max_band_count();
    if (max_bands > MAX_PRECOMP)
    {
        max_bands = MAX_PRECOMP;
    }
    if (max_bands > precomp_max_bands || old_tolerance != Pz.tol)
    {
        destroy_precomp_array();
        create_precomp_array(max_bands);

        fprintf(stderr,"precomputing sulston for %d bands, tolerance %d: start\n", max_bands,Pz.tol);
        fflush(0);

        b1_save = C1z.nbands;
        b2_save = C2z.nbands;
        m_save = Sz.match;

        for (b1 = 1; b1 <= precomp_max_bands; b1++)
        {
            C1z.nbands = b1;
            for (b2 = 1; b2 <= b1; b2++)
            {
                C2z.nbands = b2;
                for (m = 0; m <= b2; m++)
                {
                    Sz.match = m;
                    choose_sulston(&C1z, &C2z, &Sz);
                    sulston_precomp[b1][b2][m] = Sz.prob;
                }
            }
        }

        C1z.nbands = b1_save;;
        C2z.nbands = b2_save;
        Sz.match = m_save;

        fprintf(stderr,"precomputing sulston: done\n");fflush(0);
    }
    old_tolerance = Pz.tol;
}
